#include <stdexcept>
#include <sstream>

#include "../Okada.hxx"

Okada::Displacement Okada::dc3d(const FTensor::Tensor1<double,3> &coord)
{
  if(coord(2)>0)
    {
      std::stringstream ss;
      ss << "z must be less than zero, but the input was: " << coord(2);
      throw std::runtime_error(ss.str());
    }
  
  FTensor::Tensor1<double,3> u(0,0,0);
  FTensor::Tensor2<double,3,3> du(0,0,0,0,0,0,0,0,0);

  FTensor::Index<'a',3> a;
  FTensor::Index<'b',3> b;

  double xi[2]={clamp(coord(0)), clamp(coord(0)-length)};
  
  Displacement real(source_contribution(true,xi,coord));
  Displacement image(source_contribution(false,xi,coord));
  real.first(a)+=image.first(a);
  real.second(a,b)+=image.second(a,b);
  return real;
}
